A variant transfer matrix method suitable for transport through multi-probe systems 
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We have developed a variant transfer matrix method that is suitable for transport through multi- 
probe systems. Using this method, we have numerically studied the quantum spin-Hall effect 
(QSHE) on the 2D graphene with both intrinsic (V ao ) and Rashba (V r ) spin-orbit(SO) couplings. 
The integer QSHE arises in the presence of intrinsic SO interaction and is gradually destroyed by 
the Rashba SO interaction and disorder fluctuation. We have numerically determined the phase 
boundaries separating integer QSHE and spin-Hall liquid. We have found that when V ao > 0.2t 
with t is hopping constant, the energy gap needed for the integer QSHE is the largest satisfying 
\E\ < t. For smaller V ao the energy gap decreases linearly. In the presence of Rashba SO interaction 
or disorders, the energy gap diminishes. With Rashba SO interaction the integer QSHE is robust 
at the largest energy within the energy gap while at the smallest energy within the energy gap the 
integer QSHE is insensitive to the disorders. 

PACS numbers: 71.70.Ej, 72.15.Rn, 72.25.-b 
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INTRODUCTION 



Graphene is a 2-dimensional honeycomb lattice of sin- 
gle atomic carbon layer and has a special band structures. 
With more and more experimental discoveries and the- 
oretical predictions [lj, [2j, ba, |j , la, [f| , there is currently a 
intense interest on electronic properties on the graphene 
sheet. Especially the spin Hall effect (SHE) has the po- 
tential to provide a purely electrical means to control 
the spins of electron in the absence of non-ferromagnetic 
materials and magnetic field [7]. This is because the spin- 
orbit interaction in the Graphene exerts a torque on the 
spin of electron whose precessing leads to a spin polarized 
current. In a four probe device, this spin polarized cur- 
rent can lead to a pure spin current without accompany- 
ing charge current [8|. It has been proposed by Haldane[9| 
that a quantum Hall effect may exist in the absence of 
magnetic field. Similarly, integer quantum spin-Hall ef- 
fect can exist on a honeycomb lattice when the intrinsic 
spin orbit interaction is present [71 1 1 01] - In the presence 
of disorder the charge conductance of mesoscopic conduc- 
tors show universal features with a universal conductance 
fluctuation JdJ and the spin-Hall conductance also fluc- 
tuates with a universal value [12j in the presence of spin 
orbit interaction. The presence of disorder can also de- 
stroy the integer quantum spin-Hall effect and quantum 
Hall effect f 3] for a Graphene system with intrinsic spin 
orbit interaction^?]. Hence it is of interest to map out 
the phase diagram for the integer quantum spin-Hall ef- 
fect. In this paper, we investigate the disorder effect on 
the spin Hall current for a four-probe Graphene system 
in the presence of intrinsic and/or Rashba SO interac- 
tions, denoted as V so and V r , respectively. For such a 
system, the conventional transfer matrix method can not 
be used. So the direct matrix inversion method must 
be used to obtain the Green's function that is needed 
for the transport properties. As a result, the simulation 



of a multi-probe system using the direct method is very 
calculational demanding. 

In this paper, we developed an algorithm based on the 
idea of transfer matrix that is much faster than the direct 
method. As an application, we have numerically mapped 
out the phase diagram for a two dimensional honeycomb 
lattice in the presence of the intrinsic and/or Rashba SO 
interactions and disorders. When turning on the Rashba 
SO interaction, we found that the energy gap needed 
for the IQSHE is \E\ < t for V so > 0.2t and decreases 
linearly when V so < 0.2t. In the presence of Rashba 
SO interaction, the phase diagram (E, V r ) is asymmetric 
about the Fermi energy. The IQSHE is more difficult 
to destroy at the largest energy of the energy gap. In 
the presence of disorder, the phase diagram (E, W) is 
again asymmetric about the Fermi energy but it is the 
smallest energy of the energy gap that is robust against 
the disorder fluctuation. 



THEORETICAL FORMALISM 

In the tight-binding representation, the Hamiltonian 
for the 2D honeycomb lattice of the graphene structure 
can be written as[7J, [9(: 

% 



+ iV r ^2 c|e z -(crxdij)Cj+53eictci (I) 

where c](ci) is electron creation (annihilation) opera- 
tor and a are Pauli matrices. The first term is due to 
the nearest hopping. The second term is the intrinsic 
spin-orbit interaction that involves the next nearest sites. 
Here i and j are two next nearest neighbor sites, k is the 
common nearest neighbor of i and j, and djfc describes a 



vector pointing from k to i. The third term is due to the 
Rashba spin-orbit coupling. The last term is the on-site 
energy where e^ is a random on-site potential uniformly 
distributed in the interval [—W/2,W/2]. In this Hamil- 
tonian, we have set the lattice constant to be unity. 

We consider a four-probe device as shown schemati- 
cally in FIG. la. The four probes are exactly the exten- 
sion from the central scattering region, i.e., the probes 
are graphene ribbons. The number of sites in the scat- 
tering region is denoted as N = n x xn y , where there are 
n x = 8xn + 1 sites on n y = 4xn chains (FIG. la shows 
the cell for n = 1)[14|. We apply external bias volt- 
ages Vi with i = 1, 2, 3, 4 at the four different probes as 
Vi = (v/2, 0, —v/2, 0). In the presence of Rashba SO in- 
teraction, the spin is not a good quantum number. As a 
result, the spin current is not conserved using the conven- 
tional definition. Hence we switch off the Rashba SO in- 
teraction in the 2nd and 4th probes. Similar to the setup 
of Ref. 7] our setup can generate integer quantum spin 
Hall effect. The difference between the setup of Ref. 7] 
and ours is that the lead in Ref.[7|] is a square lattice 
without SO interactions while our lead is still honeycomb 
lattice with SO interactions except that the Rashba SO 
interaction has been switched off in lead 2 and 4. The use 
of the square lattice as a lead has two consequences. It 
provides additional interfacial scattering between scat- 
tering region and the lead due to the lattice mismatch 
and the mismatch in SO interactions. In addition, the 
dimension of the self-energy matrix for the square lattice 
lead with SO interaction is much smaller. The spin-Hall 
conductance G s h can be calculated from the multi-probe 
Landauer-Buttiker formula [a, 1 12J: 

G sH = (e/87r)[(T 2Tjl - T 2L1 ) - (T 2T , 3 - T 2i>3 )] (2) 

where the transmission coefficient is given by T 2ct .i = 
Tr(T 2a G r T 1 G a ) with G r ' a being the retarded and ad- 
vanced Green functions of the central disordered region 
which can be evaluated numerically. The quantities T^ 
are the linewidth functions describing coupling of the 
probes and the scattering region and are obtained by 
calculating self-energies S r due to the semi-infinite leads 
using a transfer matrices method 151]. In the following, 
our numerical data are mainly on a system with n = 8 or 
32x65 sites in the system. To fix units, throughout this 
paper, we define the Fermi-energy E, disorder strength 
W, intrinsic spin-orbit coupling V so and Rashba spin- 
orbit coupling V r in terms of the hopping energy t. 

For the four-probe device, the conventional transfer 
matrix that is suitable for two-probe devices can no 
longer be used. Below, we provide a modified transfer 
matrix method for the four-probe device. Note that the 
self-energy E r is a matrix with non-zero elements at those 
positions corresponding to the interface sites between a 
lead and the scattering region [l6|. Because evaluating 
the Green's function G r corresponds to the inversion of 
a matrix, a reasonable numbering scheme to the lattice 



sites can minimize the bandwidth of the matrix and thus 
reduce the cost of numerical computation. For example, 
to obtain the narrowest bandwidth for our system we 
partition the system into layers shown in FIG. lb so that 
there is no coupling between the next nearest layers. We 
then label each site layer by layer from the center of the 
system (see FIG. la). As a result, the matrix E — H — S r 
becomes a block tri-diagonal matrix: 



E - H 



( A x C x . 

-B 2 A 2 C 2 



V 



\ 



An-1 L-m-1 



where A n is a (128n — 56) x (128n — 56) matrix, C n 
is a (128n — 56) x (128n + 72) matrix, and B n is a 
(128n— 56) x(128n— 184) matrix. Here n = 1 corresponds 
to the innermost layer and n = m is for the outermost 
layer. A direct inversion of this block tri-diagonal matrix 
is already faster than the other labeling schemes. How- 
ever, if we are interested in the transmission coefficient, 
it is not necessary to invert the whole matrix. This is 
because the self-energies of the leads are coupled only 
to A m of the outermost layers, from Landauer-Buttiker 's 
formula it is enough to calculate the Green's function 
G r mm which satisfys the following equation, 
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where I m is a unit matrix of dimension m. In general, 
the solution Xi of the following equation with block tri- 
diagonal matrix can be easily obtained. 
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From the first row 

A\ X\ + C\ X 2 = Ri , 
we have 

X 1 +Ai 1 C 1 X 2 = A 1 - 1 R 1 . 
From the 2 nd row, 

B 2 X\ + A 2 X 2 + C 2 X% = R 2 , 



eliminating X±, we have 

(A 2 - B 2 A- 1 C 1 )X 2 + C 2 X 3 = R 2 - B 2 A~ 1 R 1 . 
This equation can be written as 



F 2 X 2 + C 2 X$ — D 2 , 



where 



F 2 = A 2 - B 2 A^ 1 C 1 ,D 2 = R 2 - B 2 A^ 1 Rl 

From the 3 rd row, 

BzX 2 + A3X3 + C3X4 = R3, 

eliminating X 2 , we have 

F3X3 + C3X4 = D3, 

where 

F 3 = A3 - B 3 F 2 - 1 C 2 ,D 3 =R 3 - B3F 2 T l D 2 . 

Therefore, we have the following recursion relation, 

F\ = A\, initial 

Fi= Ai-BtF^Xd-!, t = 2,3,-.-,m 

D\ = i?i, initial 

Di= fli-Bif^A-i, t = 2,3,---,ro 

Finally, we have 
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From the last row, we can solve for X„ 



^m — ^ m -Urn- 



We can cancel X m in the last but one equation 

X m -i = F m _ 1 (D m _\ — C m -\X m ). 

In our case, Xi = G r im and Ri — Si m I m and we are 
only interest in the solution G mm . Hence we have the 
solution 



G r = F' 

^ mm ri 



where 



Fi 



Ai, 



F t = Ai - BiFr^d-u i = 2,3,-..,m 



To test the speed of this algorithm, we have calculated 
the spin Hall conductance for the four-probe graphene 
system with different system size labeled by n on a mat- 
lab platform. The calculation is done at a fixed energy 
and for 1000 random configurations. The cpu times are 
listed in Table 1 where speed of direct matrix inversion 
and the algorithm just described are compared. We see 
that the speed up factor increases as the system size in- 
creases. For instance, for n — 8 which corresponds to 
2080 sites (amounts to a 4016 x 4016 matrix) in the scat- 
tering region, a factor of 100 is gained in speed. We note 
that in the presence of intrinsic SO interaction the cou- 
pling involves next nearest neighbor interaction. This 
is the major factor that slows down our algorithm. As 
shown in TABLE 1, for a square lattice without intrin- 
sic SO interaction but with Rashba SO interaction, the 
speed up factor is around 200 for a 40 x 40 system (matrix 
dimension 3200). The new algorithm is particular useful 
when the large number of disorder samples and different 
sample sizes are needed for the calculation of the con- 
ductance fluctuation and its scaling with size. Finally, 
we wish to mention that this algorithm also applies to 
multi-probe systems such as six-probe systems. 



NUMERICAL RESULTS 

It has been shown that in the presence of disorder or 
Rashba SO interaction the QSHE may be destroyed [7|. 
As an application of our algorithm, we study the phase 
phase boundary between regimes of the integer QSHE 
regime and the QSH liquid in the presence of disorder. 
For this purpose, we set a criteria for the QSH, i.e., if 
Gs#>0.999 we say it reaches an integer quantum spin 
Hall plateau (IQSH). Since the integer QSHE is due to 
the presence of intrinsic SOI, we first study the phase dia- 
gram of a clean sample in the absence of Rashba SOI, i.e., 
the two-component Haldane's model[9j. For this model, 
there is an energy gap within which the IQSH effect ex- 
ists. FIG. 2 depicts the phase diagram in (E,V so ) plane 
with a curve separates the integer QSHE and SHE liq- 
uid. We see that the phase diagram is symmetric about 
the Fermi energy E and the integer QSHE exists only 
for energy E < 1 that corresponds to the energy gap. 
FIG. 2 shows that the energy gap depends on the strength 
of intrinsic SO interaction. When V so > 0.2 the en- 
ergy gap is the largest between E = [—1,1] while for 
V so < 0.2, the energy gap gradually diminishes to zero 
in a linear fashion. Our numerical data show that for 
V so < 0.025 the IQSHE disappears (see FIG.2). Between 
V so = [0.025,0.18] the phase boundary is a linear curve. 
When Vso > 0.20, the phase boundary becomes a sharp 
vertical line. 

For Haldane's model, the a z is a good quantum num- 
ber. However, in the presence of Rashba SOI the spin 
experiences a spin torque while traversing the system. 



This can destroy the IQSHE at large enough Rashba 
SOI strength V r . In FIG. 3, we show the spin-Hall con- 
ductance G s h vs Fermi energy at difference V r when 
V so = 0.1,0.2. In FIG. 3a we see that when V r = 0, the 
spin-Hall conductance is quantized between E = —0.52 
and +0.52. As V r increases to 0.1, and the energy gap de- 
creases to —0.22 and 0.51. Upon further increasing V r to 
0.2 and 0.3, the gaps shrink to, respectively, [0.06,0.50] 
and [0.34,0.46]. In Ref.0] the IQSHE is completely de- 
stroyed when V r = 0.3 which is different from our re- 
sult. The difference is due to the lead used in Ref. 7] 
that causes additional scattering. The larger intrinsic 
SO interaction strength V so , the more difficult to destroy 
the integer QSHE as can be seen from FIG. 3b. 

In the presence of Rashba SO interaction the phase dia- 
gram in (E, V r ) plane at different intrinsic SO interaction 
strengths is shown in FIG. 4. We see that the phase dia- 
gram is asymmetric about the Fermi energy and it is more 
difficult to destroy the integer QSHE for largest positive 
energies within the energy gap, e.g., near E — 0.51 when 
V so = 0.1. Similar to FIG. 2, we see that when V so > 0.2 
integer QSHE can exist for all energies as long as \E\ < 1. 
Roughly speaking, the energy gap decreases linearly with 
increasing of Rashba SOI and there is a threshold V r be- 
yond which the integer QSHE disappears. For instance, 
when V r > 0.3 and V so = 0.1, the integer QSHE is de- 
stroyed. 

From the above analysis, we see that V so = 0.2 is 
an important point separating two different behaviors 
in (E,V so ) and (E,V r ) phase diagrams. Now we ex- 
amine the effect of disorder on the QSHE. FIG. 5 shows 
the phase diagram of integer QSHE on (E, W) at two 
typical intrinsic SO interaction strengths V so = 0.1 and 
V so = 0.2. The phase diagrams are asymmetric about the 
Fermi energy. Generally speaking, the larger the Rashba 
SO interaction strength V r , the smaller the energy gap 
needed for integer QSHE. We already see from FIG. 4 
that the integer QSHE is more robust against Rashba SO 
interaction strength V r at positive Fermi energy within 
the energy gap. In contrast, it is small Fermi energies 
within the energy gap that are stable against the disor- 
der fluctuation, especially for large Rashba SO interac- 
tion strength. In addition, the phase boundary at posi- 
tive Fermi energy are not very sensitive to the variation 
of Rashba SO interaction strength. The larger the in- 
trinsic SO interaction, the larger the disorder strength 
W c needed to destroy the integer QSHE. In FIG. 6, we 
estimate this critical disorder strength W c and plot it vs 
V so for E = 0.01 and V r = 0. 

If we replace the Rashba SO interaction by the Dres- 
selhaus SO interaction, we have numerically confirmed 
that the phase diagram of IQSHC in (E, W) plane is the 
same if we change E by —E. 



In summary, we have developed variant transfer matrix 
method that is suitable for multi-probe systems. With 
this algorithm, the speed gained is of a factor 100 for a 
system of 2080 sites with the next nearest SO interac- 
tion on a honeycomb lattice. For the square lattice with 
Rashba SO interaction, the speed gained is around 200 
for a 40 x 40 system. Using this algorithm, we have stud- 
ied the phase diagrams of the graphene with intrinsic and 
Rashba SO interaction in the presence of disorder. 
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FIG. 1: (Color online) Schematic plot of the four terminal 
mesoscopic sample where the intrinsic SO interaction exists 
in the center scattering region and the leads 1, 3. And the 
Rashba SO only exists in the center part and the leads 1, 3, 
when the spin-Hall conductance is measured through leads 2, 
4. 
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FIG. 2: (Color online)Phase diagram of IQSHC on (E, V ao ) 
plane for W — and V r = 0. The curve separates the IQSHC 
regime and the spin-Hall liquid regime. 
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FIG. 3: (Color online)spin Hall Conductance versus electron 
Fermi energy for V r =0, 0.1, 0.2, 0.3 on the N = 32x65 sample. 
(a) for W = and V s „ = 0.1; (b) for W = and V s0 = 0.2. 
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FIG. 4: (Color online)Phase diagram for integer quantum 
spin Hall conductance on (E,V r ) plane. Squares, circles, left- 
triangles and right-triangles are for l/ so =0.1, 0.2, 0.3, 0.4, re- 
spectively. The areas encircled by the curves and the V r —0 
line are the integer quantum spin Hall conductance regimes 
for different intrinsic SOI. 
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FIG. 5: (Color online) Phase diagram of IQSHC on (E, W) 
plane for different Rashba SO coupling in the presence of (a) 
V so = 0.1 and (b)V ao = 0.2. Squares, circles, stars and nimbus 
are for V r =0, 0.1, 0.2, 0.3. The areas encircled by the curves 
and the W—0 line are the IQSHC regimes for different Rashba 
SOI. 
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FIG. 6: The critical disorder strength versus intrinsic SO cou- 
pling Vso- The corresponding Fermi energy is E — 0.01 and 
V r — 0. The spin Hall conductance in the regime encircled by 
the curve and the Rashba SOI axis is well quantized. 



Tablel 



a four-probe graphene 


a four-probe square lattice 


system size n 


direct method 


our method 


system size L 


direct method 


our method 


6 


14803s 


383s 


30 


5102s 


53s 


7 


31491s 


606s 


40 


29337s 


152s 


8 


104275s 


958s 


50 


89917s 


309s 



TABLE I: The cpu times for different system sizes using dif- 
ferent methods are calculated at a fixed Fermi energy for 1000 
random configurations. 



